function M = regiongr2( I, x, y)

[SizeY, SizeX] = size(I);

%konstanty
ka = 0.05;
kb = 0.6;
kx = 20 / (SizeX * SizeY);  
ky = 0.6*kx;

%inicializace
M = zeros(SizeY, SizeX);  %oblast

xi = floor(x); 
xr = x - xi;
yi = floor(y);
yr = y - yi;
M(yi, xi) = (1-xr) * (1-yr); %pocatek
M(yi, xi + 1) = (1-xr) * yr;
M(yi + 1, xi) = (1-yr) * xr;
M(yi + 1, xi + 1) = xr * yr;

Dx = kx * (((1:SizeX) - x) .^ 2);
Dy = ky * (((1:SizeY) - y) .^ 2);
D1 = repmat(Dx,SizeY,1) + repmat(Dy',1,SizeX);   %vzdalenosti od stredu
Dx = kx * (((1:SizeX) - SizeX/2) .^ 2);
Dy = ky * (((1:SizeY) - SizeY/2) .^ 2);
D2 = repmat(Dx,SizeY,1) + repmat(Dy',1,SizeX);
D = 1 + D1 + 0.5*D2;

MsizePrev = 0;    %velikost oblasti v predchozi iteraci
Msize = 1;        %nova velikost oblasti
iter = 0;         %pocet iteraci

while (Msize - MsizePrev > ka * iter)        
    iter = iter + 1;
    mean = sum(sum(M .* I)) / Msize;  %prumerna barva oblasti
      
    R = abs(I - mean);     %rozdily mezi barvou pixelu a barvou oblasti
    
    Q = kb - R .* D;  %koeficienty pixelu
    Q(Q < 0) = 0;
    
    S = zeros(SizeY, SizeX);   %soucety sousednich hodnot
    S(1:(SizeY - 1),:) = M(2:SizeY,:);
    S(2:SizeY,:) = S(2:SizeY,:) + M(1:(SizeY - 1),:);
    S(:,1:(SizeX - 1)) = S(:,1:(SizeX - 1)) + M(:, 2:SizeX);
    S(:,2:SizeX) = S(:,2:SizeX) + M(:,1:(SizeX - 1));
    
    M = M + S .* Q;    %zvetseni oblasti
    M(M > 1) = 1;
    
    MsizePrev = Msize;
    Msize = sum(sum(M)); 
end
M = (M == 1); 
end

